Assignment 10

In [1]:

    

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from time import time
%matplotlib inline


    

In [2]:

    

def Runge_Kutta(f, t, x0):
    n = len(x0)
    N = len(t)
    x = np.zeros((n, N))
    x[:, 0] = x0
    dt = t[1] - t[0]
    for i in range(N - 1):
        k1 = f(x[:, i])
        k2 = f(x[:, i] + dt * k1 / 2)
        k3 = f(x[:, i] + dt * k2 / 2)
        k4 = f(x[:, i] + dt * k3)
        x[:, i + 1] =  x[:, i] + dt * (k1 + 2 *  (k2 + k3) + k4) / 6
    return x


    

In [3]:

    

def Lorenz(X, t=0):
    return np.array([10 * (X[1] - X[0]),
                     28 * X[0] - X[1] - X[0] * X[2],
                     X[0] * X[1] - 8 * X[2] / 3])


    

In [4]:

    

N = 10 ** 4
t = np.linspace(0, 50, N)
x0 = np.array([1, 0, 0])
x = Runge_Kutta(Lorenz, t, x0)
plt.plot(t, x.T)
plt.xlabel('t')
plt.ylabel('x, y, z')


    

    
Out[4]:





<matplotlib.text.Text at 0x7f59cc4eeb38>






    

In [5]:

    

def Lorenz100(X, t=0):
    return np.hstack((10 * (X[100:200] - X[:100]),
                      28 * X[:100] - X[100:200] - X[:100] * X[200:300],
                      X[:100] * X[100:200] - 8 * X[200:300] / 3))


    

In [6]:

    

N = 10 ** 5
t = np.linspace(0, 50, N)
x0 = np.hstack((np.random.normal(1, .05, size=100), np.zeros(200)))
current_time = time()
x = Runge_Kutta(Lorenz100, t, x0)
print('It took', time() - current_time, 's')


    

    




It took 71.88326144218445 s

In [7]:

    

x1 = np.mean(x[:100, :], axis=0)
x2 = np.mean(x[100:200, :], axis=0)
x3 = np.mean(x[200:300, :], axis=0)


    

In [8]:

    

plt.plot(t, x1)
plt.plot(t, x2)
plt.plot(t, x3)
plt.xlabel('t')
plt.ylabel('mean of x, mean of y, mean of z')


    

    
Out[8]:





<matplotlib.text.Text at 0x7f598f0a0c50>






    

In [9]:

    

var_x1 = np.var(x[:100, :], axis=0)
var_x2 = np.var(x[100:200, :], axis=0)
var_x3 = np.var(x[200:300, :], axis=0)

plt.plot(t, var_x1)
plt.plot(t, var_x2)
plt.plot(t, var_x3)
plt.ylabel('var of x, var of y, var of z')


    

    
Out[9]:





<matplotlib.text.Text at 0x7f598f0a5240>






    

In [10]:

    

def J(x,y):
    return np.array([[-.1 * y, -.1 * x + 3 * y], [-.5 * y, -.5 * x]])

print('The eigenvalues in sqrt(15), 1 / sqrt(15):')
print(np.linalg.eigvals(J(np.sqrt(15), 1 / np.sqrt(15))))
print('\nThe eigenvalues in sqrt(15), 1 / sqrt(15):')
print(np.linalg.eigvals(J(-np.sqrt(15), -1 / np.sqrt(15))))


    

    




The eigenvalues in sqrt(15), 1 / sqrt(15):
[-0.05235727 -1.90995429]

The eigenvalues in sqrt(15), 1 / sqrt(15):
[ 0.05235727  1.90995429]

In [11]:

    

J = np.array([[-10, 10, 0], [28, -1, 0], [0, 0, -8 / 3]])
print('The eigenvalues in 0, 0, 0:')
print(np.linalg.eigvals(J))


    

    




The eigenvalues in 0, 0, 0:
[-22.82772345  11.82772345  -2.66666667]

In [ ]: